import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from scipy import constants
a0 = constants.physical_constants['Bohr radius'][0]
Angstrom = constants.angstrom
Here is a function to calculate the spherical coordinates r , ϕ , θ of a point with Cartesian coordinates x , y , z .
def r_phi_theta(x, y, z):
r = np.sqrt(x**2 + y**2 + z**2)
phi = np.arctan2(y, x) + np.pi
theta = np.arccos(z/r)
return r, phi, theta
The next cell uses numpy.mgrid to create a 3d grid x y z grid spanning [ − 20 , 20 ] Angstroms in each dimension. The grid has 48
110592 gridpoints. More gridpoints means better-looking graphics but also longer compute times. mgrid gives us three individual arrays: one contains the x coordinates of all points in our 3d grid, another has the y coordinates, and the third has the z coordinates.
The second line of code uses the r_phi_theta function defined above to calculate three arrays that give us the same 3d grid as above, but in spherical coordinates.
You will notice how the output values of both mgrid and r_phi_theta() are "unpacked" and assigned to three separate values. Non-unpacked syntax such as
grid = r_phi_theta(X, Y, Z) would return a tuple containing the three spherical coordinate grids.
X, Y, Z = np.mgrid[-20*a0:20*a0:48j, -20*a0:20*a0:48j, -20*a0:20*a0:48j]
r_grid, phi_grid, theta_grid = r_phi_theta(X, Y, Z)
The next cell contains a probability function that performs the operation r 2 ψ ∗ ( r , ϕ , θ ) ψ ( r , ϕ , θ ) . The second line normalizes the probability function in a way that makes it easy to plot. The requirement for the probability function is that ∫ ∞
0 ∫ 2 π
0 ∫ π
0 r 2 ψ ( r , ϕ , θ ) ∗ ψ ( r , ϕ , θ ) d r d ϕ d
1. However, since the volume of the hydrogen atom is so tiny (recall that the Bohr radius is less than an Angstrom), the probability densities can be huge, of order 10 26 m − 3 . We will put our probability densities on a scale of 0 to 1 just so the numbers on our colorbars won't look insane.
def probfunc(r, wavefunc):
p = r**2 * np.conj(wavefunc) * wavefunc
return p.real / np.max(p.real)
# The .real syntax picks off the real part of a complex number. Even though the
# imaginary part of the probability is zero, variable p is still represented in
# the computer as a complex number. Our plotting routines will only work on
# arrays of real numbers.
The next cell defines two wave functions: n1l0m0(r) for n , l ,
1 , 0 , 0 and n2l0m0(r) for n , l ,
2 , 0 , 0 .
How come these wave functions do not depend on ϕ and θ ?
n1l0m0 = lambda r: np.exp(-r/a0) / (np.sqrt(np.pi) * a0**1.5)
n2l0m0 = lambda r: (2 - r/a0) * np.exp(-r/(2*a0)) / (4 * np.sqrt(2*np.pi) * a0**1.5)
prob_n1l0m0 = probfunc(r_grid, n1l0m0(r_grid))
prob_n2l0m0 = probfunc(r_grid, n2l0m0(r_grid))
Now we use plotly to plot isosurfaces of our probability functions An isosurface is like a contour line on a map, except in 3d: it's a surface over which the value of the dependent variable is constant. These visualizations will each have eight semi-transparent isosurfaces. Darker surfaces correspond to higher probability densities. You can click and drag to rotate, enlarge, or shrink your visualization.
You will see from the graphics that n , l ,
1 , 0 , 0 and n , l ,
2 , 0 , 0 give spherically symmetric probability functions, as expected.
n =1 ,
0 , m
0
# n = 1, l = 0, m = 0
fig = go.Figure(data=go.Isosurface(
x=X.flatten() / Angstrom, # x values of the grid in Angstroms
y=Y.flatten() / Angstrom, # y values of the grid in Angstroms
z=Z.flatten() / Angstrom, # z values of the grid in Angstroms
value=prob_n1l0m0.flatten(), # independent variable
isomin=0.05, # Minimum normalized probability density to render in an isosurface
isomax=0.95, # Maximum normalized probability density to render in an isosurface
opacity=0.4, # Set a low opacity so each surface is partially transparent
colorscale='Plotly3_r', # Nice-looking color table
surface_count=8, # number of isosurfaces to plot (2 by default: only min and max)
colorbar_nticks=8, # colorbar ticks correspond to isosurface values
caps=dict(x_show=False, y_show=False, z_show=False))
)
# Change axis lables and make the plot larger than the default
fig.update_layout(scene = dict(
xaxis_title='x (Angstroms)',
yaxis_title='y (Angstroms)',
zaxis_title='z (Angstroms)'),
width=700,
margin=dict(r=10, b=10, l=10, t=10))
fig.show()
# n = 2, l = 0, m = 0
fig = go.Figure(data=go.Isosurface(
x=X.flatten() / Angstrom, # x values of the grid in Angstroms
y=Y.flatten() / Angstrom, # y values of the grid in Angstroms
z=Z.flatten() / Angstrom, # z values of the grid in Angstroms
value=prob_n2l0m0.flatten(), # independent variable
isomin=0.05, # Minimum normalized probability density to render in an isosurface
isomax=0.95, # Maximum normalized probability density to render in an isosurface
opacity=0.4, # Set a low opacity so each surface is partially transparent
colorscale='Plotly3_r', # Nice-looking color table
surface_count=8, # number of isosurfaces to plot (2 by default: only min and max)
colorbar_nticks=8, # colorbar ticks correspond to isosurface values
caps=dict(x_show=False, y_show=False, z_show=False)
))
# Change axis lables and make the plot larger than the default
fig.update_layout(scene = dict(
xaxis_title='x (Angstroms)',
yaxis_title='y (Angstroms)',
zaxis_title='z (Angstroms)'),
width=700,
margin=dict(r=10, b=10, l=10, t=10))
fig.show()
Our isosurface renderings might not help us visualize a probability density that is not monotonically increasing as a function of radius: a low probability region (light pink) might be concealed insided a higher-probability region (purple). To visualize the radial structure of our probability functions in more detail, we add a slice to each isosurface rendering. The slice I chose is parallel to the y − z plane and centered at
0 . See if you can spot the one-line change to the code that adds the slice.
n =1 ,
0 ,
0
# n = 1, l = 0, m = 0
fig = go.Figure(data=go.Isosurface(
x=X.flatten() / Angstrom, # x values of the grid in Angstroms
y=Y.flatten() / Angstrom, # y values of the grid in Angstroms
z=Z.flatten() / Angstrom, # z values of the grid in Angstroms
value=prob_n1l0m0.flatten(), # independent variable
isomin=0.05, # Minimum normalized probability density to render in an isosurface
isomax=0.95, # Maximum normalized probability density to render in an isosurface
opacity=0.4, # Set a low opacity so each surface is partially transparent
colorscale='Plotly3_r', # Nice-looking color table
surface_count=8, # number of isosurfaces to plot (2 by default: only min and max)
colorbar_nticks=8, # colorbar ticks correspond to isosurface values
slices_x=dict(show=True, locations=[0]),
caps=dict(x_show=False, y_show=False, z_show=False))
)
# Change axis lables and make the plot larger than the default
fig.update_layout(scene = dict(
xaxis_title='x (Angstroms)',
yaxis_title='y (Angstroms)',
zaxis_title='z (Angstroms)'),
width=700,
margin=dict(r=10, b=10, l=10, t=10))
fig.show()
# n = 2, l = 0, m = 0
fig = go.Figure(data=go.Isosurface(
x=X.flatten() / Angstrom, # x values of the grid in Angstroms
y=Y.flatten() / Angstrom, # y values of the grid in Angstroms
z=Z.flatten() / Angstrom, # z values of the grid in Angstroms
value=prob_n2l0m0.flatten(), # independent variable
isomin=0.05, # Minimum normalized probability density to render in an isosurface
isomax=0.95, # Maximum normalized probability density to render in an isosurface
opacity=0.4, # Set a low opacity so each surface is partially transparent
colorscale='Plotly3_r', # Nice-looking color table
surface_count=8, # number of isosurfaces to plot (2 by default: only min and max)
colorbar_nticks=8, # colorbar ticks correspond to isosurface values
caps=dict(x_show=False, y_show=False, z_show=False),
slices_x=dict(show=True, locations=[0]),
))
# Change axis lables and make the plot larger than the default
fig.update_layout(scene = dict(
xaxis_title='x (Angstroms)',
yaxis_title='y (Angstroms)',
zaxis_title='z (Angstroms)'),
width=700,
margin=dict(r=10, b=10, l=10, t=10))
fig.show()
Based on what you see in the slices, how are the probability functions for
1 ,
0 ,
0 and
2 ,
0 ,
0 different?
Try
2 ,
1 ,
0
n2l1m0 = lambda r, theta: (r/a0) * np.exp(-r/(2*a0)) * np.cos(theta) \
/ (4 * np.sqrt(2*np.pi) * a0**1.5)
# n = 2, l = 1, m = 0
prob_n2l1m0 = probfunc(r_grid, n2l1m0(r_grid, theta_grid))
fig = go.Figure(data=go.Isosurface(
x=X.flatten() / Angstrom, # x values of the grid in Angstroms
y=Y.flatten() / Angstrom, # y values of the grid in Angstroms
z=Z.flatten() / Angstrom, # z values of the grid in Angstroms
value=prob_n2l1m0.flatten(), # independent variable
isomin=0.05, # Minimum normalized probability density to render in an isosurface
isomax=0.95, # Maximum normalized probability density to render in an isosurface
opacity=0.4, # Set a low opacity so each surface is partially transparent
colorscale='Plotly3_r', # Nice-looking color table
surface_count=8, # number of isosurfaces to plot (2 by default: only min and max)
colorbar_nticks=8, # colorbar ticks correspond to isosurface values
caps=dict(x_show=False, y_show=False, z_show=False),
# slices_x=dict(show=True, locations=[0]),
))
# Change axis lables and make the plot larger than the default
fig.update_layout(scene = dict(
xaxis_title='x (Angstroms)',
yaxis_title='y (Angstroms)',
zaxis_title='z (Angstroms)'),
width=700,
margin=dict(r=10, b=10, l=10, t=10))
fig.show()
n2l1m1 = lambda r, phi, theta: (r/a0) * np.exp(-r/(2*a0)) * np.sin(theta) \
* np.exp(1j * phi) / (8 * np.sqrt(3*np.pi) * a0**1.5)
# n = 2, l = 1, m = 1
prob_n2l1m1 = probfunc(r_grid, n2l1m1(r_grid, phi_grid, theta_grid))
fig = go.Figure(data=go.Isosurface(
x=X.flatten() / Angstrom, # x values of the grid in Angstroms
y=Y.flatten() / Angstrom, # y values of the grid in Angstroms
z=Z.flatten() / Angstrom, # z values of the grid in Angstroms
value=prob_n2l1m1.flatten(), # independent variable
isomin=0.05, # Minimum normalized probability density to render in an isosurface
isomax=0.95, # Maximum normalized probability density to render in an isosurface
opacity=0.4, # Set a low opacity so each surface is partially transparent
colorscale='Plotly3_r', # Nice-looking color table
surface_count=8, # number of isosurfaces to plot (2 by default: only min and max)
colorbar_nticks=8, # colorbar ticks correspond to isosurface values
caps=dict(x_show=False, y_show=False, z_show=False),
# slices_x=dict(show=True, locations=[0]),
))
# Change axis lables and make the plot larger than the default
fig.update_layout(scene = dict(
xaxis_title='x (Angstroms)',
yaxis_title='y (Angstroms)',
zaxis_title='z (Angstroms)'),
width=700,
margin=dict(r=10, b=10, l=10, t=10))
fig.show()
# n = 2, l = 1, m = 1
fig = go.Figure(data=go.Isosurface(
x=X.flatten() / Angstrom, # x values of the grid in Angstroms
y=Y.flatten() / Angstrom, # y values of the grid in Angstroms
z=Z.flatten() / Angstrom, # z values of the grid in Angstroms
value=prob_n2l1m1.flatten(), # independent variable
isomin=0.05, # Minimum normalized probability density to render in an isosurface
isomax=0.95, # Maximum normalized probability density to render in an isosurface
opacity=0.4, # Set a low opacity so each surface is partially transparent
colorscale='Plotly3_r', # Nice-looking color table
surface_count=8, # number of isosurfaces to plot (2 by default: only min and max)
colorbar_nticks=8, # colorbar ticks correspond to isosurface values
caps=dict(x_show=False, y_show=False, z_show=False),
slices_x=dict(show=True, locations=[0]),
))
# Change axis lables and make the plot larger than the default
fig.update_layout(scene = dict(
xaxis_title='x (Angstroms)',
yaxis_title='y (Angstroms)',
zaxis_title='z (Angstroms)'),
width=700,
margin=dict(r=10, b=10, l=10, t=10))
fig.show()
# n = 2, l = 1, m = 1
fig = go.Figure(data=go.Isosurface(
x=X.flatten() / Angstrom, # x values of the grid in Angstroms
y=Y.flatten() / Angstrom, # y values of the grid in Angstroms
z=Z.flatten() / Angstrom, # z values of the grid in Angstroms
value=prob_n2l1m1.flatten(), # independent variable
isomin=0.05, # Minimum normalized probability density to render in an isosurface
isomax=0.95, # Maximum normalized probability density to render in an isosurface
opacity=0.4, # Set a low opacity so each surface is partially transparent
colorscale='Plotly3_r', # Nice-looking color table
surface_count=8, # number of isosurfaces to plot (2 by default: only min and max)
colorbar_nticks=8, # colorbar ticks correspond to isosurface values
caps=dict(x_show=False, y_show=False, z_show=False),
slices_z=dict(show=True, locations=[0]),
))
# Change axis lables and make the plot larger than the default
fig.update_layout(scene = dict(
xaxis_title='x (Angstroms)',
yaxis_title='y (Angstroms)',
zaxis_title='z (Angstroms)'),
width=700,
margin=dict(r=10, b=10, l=10, t=10))
fig.show()
n3l2m1 = lambda r, phi, theta: (r**2/a0**2) * np.exp(-r/(3*a0)) * np.sin(theta) * \
np.cos(theta) * (np.cos(phi) + 1j * np.sin(phi)) / \
(81 * np.sqrt(np.pi) * a0**1.5)
# n = 3, l = 2, m = 1
prob_n3l2m1 = probfunc(r_grid, n3l2m1(r_grid, phi_grid, theta_grid))
fig = go.Figure(data=go.Isosurface(
x=X.flatten() / Angstrom, # x values of the grid in Angstroms
y=Y.flatten() / Angstrom, # y values of the grid in Angstroms
z=Z.flatten() / Angstrom, # z values of the grid in Angstroms
value=prob_n3l2m1.flatten(), # independent variable
isomin=0.05, # Minimum normalized probability density to render in an isosurface
isomax=0.95, # Maximum normalized probability density to render in an isosurface
opacity=0.4, # Set a low opacity so each surface is partially transparent
colorscale='Plotly3_r', # Nice-looking color table
surface_count=8, # number of isosurfaces to plot (2 by default: only min and max)
colorbar_nticks=8, # colorbar ticks correspond to isosurface values
caps=dict(x_show=False, y_show=False, z_show=False),
))
# Change axis lables and make the plot larger than the default
fig.update_layout(scene = dict(
xaxis_title='x (Angstroms)',
yaxis_title='y (Angstroms)',
zaxis_title='z (Angstroms)'),
width=700,
margin=dict(r=10, b=10, l=10, t=10))
fig.show()
# n = 3, l = 2, m = 1
prob_n3l2m1 = probfunc(r_grid, n3l2m1(r_grid, phi_grid, theta_grid))
fig = go.Figure(data=go.Isosurface(
x=X.flatten() / Angstrom, # x values of the grid in Angstroms
y=Y.flatten() / Angstrom, # y values of the grid in Angstroms
z=Z.flatten() / Angstrom, # z values of the grid in Angstroms
value=prob_n3l2m1.flatten(), # independent variable
isomin=0.05, # Minimum normalized probability density to render in an isosurface
isomax=0.95, # Maximum normalized probability density to render in an isosurface
opacity=0.4, # Set a low opacity so each surface is partially transparent
colorscale='Plotly3_r', # Nice-looking color table
surface_count=8, # number of isosurfaces to plot (2 by default: only min and max)
colorbar_nticks=8, # colorbar ticks correspond to isosurface values
caps=dict(x_show=False, y_show=False, z_show=False),
slices_x=dict(show=True, locations=[0]),
))
# Change axis lables and make the plot larger than the default
fig.update_layout(scene = dict(
xaxis_title='x (Angstroms)',
yaxis_title='y (Angstroms)',
zaxis_title='z (Angstroms)'),
width=700,
margin=dict(r=10, b=10, l=10, t=10))
fig.show()
You may use one of the wave functions tabulated at the top of this assignment or find a more complex one if you're interested in intricate shapes.
n =3 , ℓ =0 , m =0
n3l0m0 = lambda r: np.exp(-r/(2*a0)) * (27 - (18*r/a0) + (2*r**2/a0**2)) / (81 * np.sqrt(3 * np.pi) * a0**1.5)
# n = 3, l = 0, m = 0
prob_n3l0m0 = probfunc(r_grid, n3l0m0(r_grid))
fig = go.Figure(data=go.Isosurface(
x=X.flatten() / Angstrom, # x values of the grid in Angstroms
y=Y.flatten() / Angstrom, # y values of the grid in Angstroms
z=Z.flatten() / Angstrom, # z values of the grid in Angstroms
value=prob_n3l0m0.flatten(), # independent variable
isomin=0.05, # Minimum normalized probability density to render in an isosurface
isomax=0.95, # Maximum normalized probability density to render in an isosurface
opacity=0.4, # Set a low opacity so each surface is partially transparent
colorscale='Plotly3_r', # Nice-looking color table
surface_count=8, # number of isosurfaces to plot (2 by default: only min and max)
colorbar_nticks=8, # colorbar ticks correspond to isosurface values
caps=dict(x_show=False, y_show=False, z_show=False),
))
# Change axis lables and make the plot larger than the default
fig.update_layout(scene = dict(
xaxis_title='x (Angstroms)',
yaxis_title='y (Angstroms)',
zaxis_title='z (Angstroms)'),
width=700,
margin=dict(r=10, b=10, l=10, t=10))
fig.show()